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Coherent propagation of two interacting particles in Id weak random potential is considered. An 
accurate estimate of the matrix element of interaction in the basis of localized states leads to mapping 
onto the relevant matrix model. This mapping allows to clarify the mechanism of enhancement of the 
localization length which turns out to be rather different from the one considered in the literature. 
Although the existence of enhancement is transparent, an analytical solution of the matrix model 
was found only for very short samples. For a more realistic situation numerical simulations were 
performed. The result of these simulations is consistent with 

h/h ~ IJ , 

where h and h are the single and two particle localization lengths and the exponent 7 depends 
on the strength of the interaction. In particular, in the limit of strong particle-particle interaction 
there is no enhancement of the coherent propagation at all (£2 ~ h)- 

PACS numbers: 72.15.Rn, 71.30.+h, 05.45. +b 



I. INTRODUCTION 

The enhancement of the propagation length for two interacting particles (TIP) in a one-dimensional random poten- 
tial was predicted in a paper by D. Shepelyansky pi a couple of years ago. This result has attracted broad interest and 



stimulated active analytical [^|J^-|9|,[l^] and numerical |^-|q,[lo|] investigations. 2-3-dimensional and quasi-ld extensions 
of the model [jD, as well as many other related problems have also been considered in various papers (see e.g. P,p|,pd| ) . 
Moreover some of these new results may even be better established than the original one jjj (see also |l3) ) . 
More specifically, the first estimate (T) of the two-particle localization length was 

h~l\ , (1) 

where fa and I2 are the one and two particle localization lengths. 

Numerical simulations first of all have confirmed the existence of enhancement. As for the specific form of fa, some 
authors || see deviations from (|l|) , whereas others || report on the agreement with ([!]) (more concretely the authors 
of H agree with |l| in the fa dependence of h although they still disagree in the overall normalization constant in 
([!]))■ However, the actual enhancement h/h which was observed in numerical simulations, turns out to be rather 
small and varies from h/h ~ 2 in the first paper [jjj to h/h ~ 3 — 6 for the most advanced computations |5]. 

We are going to consider particles moving in an exactly Id, weak, random potential. This means that the Anderson 
localization length h is much larger than the de Broglie wave length A (for the Anderson hopping Hamiltonian (|5|) 
below A ~ 1 and fa^> 1). Let the total length of the sample be L 3> fa. The main achievement of was the prediction 
of the existence of very unusual bound states for TIP (we will call them coherent states in order to distinguish from the 
usual molecular states H]). The typical distance between particles for these new states is rather large, x\ — X2 ~ fa, 
but the joint propagation length for two particles turns out to be even much larger A(x± + X2) ~ fa ^ fa- The total 
number of these new eigenstates is also sufficiently large: AiV c ~ (L/X)(fa/X). 

The central point for all the methods applied to calculation of the coherent propagation length fa was an estimate 
of the matrix element of the interaction in the basis of products of single particle localized states. However, as we will 
show in the following section, the estimate of the matrix element in ref. pi] crucially depends on the oversimplified 
assumption regarding the behavior of the single particle wave function and therefore is irrelevant. Surprisingly, all 
the authors of the following papers have accepted the matrix element estimate of jlj without any critical analysis. 

This is why the main part of our paper will be devoted to an accurate estimate of the matrix element. This estimate 
allows us to perform mapping of the two-particle problem onto the physically relevant matrix model. 
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Surely the mapping itself is impossible without some (properly motivated) assumptions and order of magnitude 
estimates. Moreover, although we see a mechanism which should lead to the enhancement of the coherent propagation 
length, we have no rigorous proof that one could not find another source of enhancement. Therefore in the light of 
contradiction with existing predictions it should be very useful to find at least one new rigorous result. To this end we 
consider in Section 3 the strong coupling limit of the Shepelyansky model. It will be shown without any assumptions 
that for the very strong interaction between particles all the enhancement disappears and I2 — h- Thus the problem 
of two interacting particles exhibits some kind of duality between the parameters U/t for weak (U <C t) and t/U 
for strong (t <C U) coupling cases |15|, where U is the strength of interaction and t is the typical single-particle 
energy In both weak and strong coupling limits the enhancement disappears. As a function of U the coherent 

propagation length I2 should reach the maximum at some value of U/t ~ 1. 

As we will show in the following section, the two main features of the matrix element of interaction lead us to 
associate the TIP problem with the matrix models which differs strongly from the previously considered ones ■ 
First of all it is the hidden hierarchy of the matrix elements. In general the matrix element of interaction in the basis 
of noninteracting two-particle states turns out to be much smaller than it was expected from the first estimate of ref. 
[ 1 1 . Only very small part of the matrix elements (namely ~ 1/ y/h) may exceed the original estimate by Shepelyansky 
This will be may be the most surprising result of our paper that such a few large matrix elements still may 
lead to a considerable enhancement. Another important feature of the matrix elements is the hidden order in their 
distribution. As we will show in Sections 5,6, after the proper ordering of the noninteracting basis the complicated 
hierarchy of the matrix elements may be described by the corresponding enveloping function of the relevant banded 
Gaussian random matrix model. 

The matrix model we are going to consider turns out to be much more complicated than those investigated in 
p|,p[-p|JlC|] . Therefore at first it would be useful to consider a simplified version of the problem. To this end in Section 
5 we consider the TIP in the short sample with the total size L being of the same order of magnitude as the single 
particle localization length L ~ l\. In this case the corresponding random matrix is the so called Power-law Random 
Band Matrix (PRBM). The elements of this matrix decrease in a power-law fashion My ~ \i— as one goes farther 
from the diagonal. Among the PRBM matrices of the general form with Mij ~ \i — j\~ a those with a = 1 correspond 
to the phase transition from localized to delocalized regime pj^ . For short sample we consider the so called inverse 
participation ratio U pr (which is effectively the number of noninteracting two-particle states mixed into one chaotic 
eigenstate). The result reads 

lipr ^ l\ j (2) 

where 7i pr is a function of the strength of the interaction approaching the value ji pr = in the weak and strong 
coupling limit. Equation (^) is proved analytically by the Renormalization Group method for ji pr -C 1 but still for 
7iprhi(Zi) 1. The basic idea of calculating (||) follows the method of Levitov JlTj who considered an even more 
complicated 3<i problem. Numerical simulations also support the result (||) for ji pr ~ 1. 

Unfortunately, for the realistic model with L li we could not find any convincing analytical solution (the matrix 
model itself will be described in Section 6). Nevertheless, we at least see that for a large sample the chaotic mixing of 
noninteracting two-particle states is systematically enhanced compared to a short sample. For example, it may be the 
same expression (||) for li pr but with 7' > 7. The perturbation theory for large L may be used only if (li pr ) — 1 <C 1. 
Nevertheless the physical expansion parameter for this perturbation theory still is \n{l{)U/t for the weak coupling 
and hx{li)t/U for the strong coupling limit. 

Therefore we have to perform numerical simulations for the matrix model associated with two particles on a large 
sample. The details of numerical procedure will be considered in Section 6 and now we will only make some general 
comments. The inverse participation ratio li pr characterizes the complexity of the typical eigenstate in the Hilbert 
space formed by noninteracting two-particle states. On the other hand the coherent propagation length I2 directly 
measures the spatial (lattice sites) dimension of the corresponding wave function. For practical calculations (Section 
6) we define I2 throw the mean squared size of the wave packet in longitudinal direction. We present our final result 
in the form 

h/h~ll ■ (3) 

First of all, we definitely see that the exponent 7 is a function of the strength of interaction. Unfortunately we do 
not know too much about this dependence of 7 on U . The only claim we can make is that 7 should go to zero for 
the weak {U -C t) and strong (t <C U) coupling limiting cases (see also Jl^l). For arbitrary U ~ t we expect 7 ~ 1. 
For the concrete choice of parameters which we have considered numerically it was also 7 < 1. However, we do not 
know whether the same inequality 7 < 1 holds for any strength of interaction. Also for l\ 3> A the exponent 7 in (j^) 
does not depend strongly on the strength of disorder or, in other words, 7 does not depend on l\. The two-particle 
localization length of the form (||) with 7 < 1 was found also numerically in || by the transfer matrix method. 



2 



Of course the expression (||) (as well as @) should be used only for large l\. This means that if one tries to extend 
(||) to small l\, 7 will become a function of l\ as well. For example, in order to fit the results of simulations, we use 
7 = 70 + cjl\. Our data (see figs. 2,3 below) shows that at least for the coherent propagation length I2 the value of 
7 does not vary strongly with l\ . However, our numerical accuracy still is not enough to exclude completely some 
exotic dependence of 7 on l\ for large 1%, say 7 = + &ln(Zi) . 

This paper is organized as follows. In Section 2 we give a general formulation of the problem and make a rough 
estimate of the matrix element of the interaction. It is shown by considering the modified Thouless block picture 
|pO|,p| that effectively the interaction between particles is enhanced by a factor ~ ln(Zi). In principle the material of 
Section 2 allows one to perform the mapping onto the matrix model which is done in Sections 5,6. However, in the 
next two sections we try to build a more stable foundation for this mapping. In Section 3 we consider the strong 
coupling limit of the Shepelyansky model. This consideration provides us with a better understanding as to how 
to distinguish regular and chaotic effects due to inter-particle interaction. On the other hand, the exact solution 
of the two-particle problem in the strong coupling regime allows us to perform the critical revision of the existing 
estimates of the coherent propagation length. In Section 4 we present a more rigorous estimate of the interaction 
matrix clement than the one in Section 2. By averaging over the disordered potential in the "unimportant" part of 
its Fourier spectrum, which is not responsible for the global features of the localized single-particle states, we get 
rid of the problem of rapid oscillations in the matrix element. In Section 5 we investigate the effect of interaction 
in the short sample. Within the Renormalization Group approach of |l7j a nontrivial solution of the model is found 
analytically at least for the weak effective interaction case. Finally, in the last Section 6 we describe the mapping of 
our coherent propagation problem onto the eigenvalue problem for some special Random Band Matrices. Features of 
this last model are investigated numerically. 



II. FORMULATION OF THE PROBLEM AND PRELIMINARY ESTIMATES 

Following Jl]] consider two particles on a Id lattice with the Hamiltonian 

H tot = H(ni) + H(n 2 )+V int . (4) 

There are no serious contradictions against considering the same problem in continuous l<i-space. The Anderson lattice 
Hamiltonian, which is used traditionally, only simplifies the numerical calculations. The single particle Hamiltonian 
has the following nonzero elements for transitions between n-th and m-th sites of the lattice 

H nm = —t{Sn, m+l + ^n,m-l) + W n 5 n ,m ■ (5) 

Here w n is the random Gaussian potential and we suppose that the disorder is weak compared to the kinetic energy: 

w n w m = w 2 S nm , w <C t . (6) 

In many papers devoted to the TIP derealization problem the disordered potential is usually chosen to be uniformly 
distributed within the range —W/2 < w n < W/2 . The comparison with our results in this case may be made via 
trivial replacement w 2 = W 2 /12 . The Hubbard on-site interaction is defined as 

(V m t^) nin2 =US ni: , l2 ^(n 1 ,n 2 ) . (7) 

In general we assume that the interaction strength U is of the same order of magnitude as the hopping matrix 
element t, though the more or less clear analytical results may be obtained only in the weak (U <§; t) and strong 
(U 3> t) coupling limits. 

For simplicity we consider distinguishable particles, but all our results are equally valid for bosons or fermions with 
opposite spins. 

The Anderson localization length has the form [|TJ 

At 2 - e 2 

h = 2^-^ , (8) 

where e is the single particle energy (— 2t < e < 2t). We will not be interested in the edges of the energy zone. 
Therefore for our purposes it will be enough to remember that l\ ss (t/w) 2 . For the lattice model both and 
li are naturally dimensionless. 
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An important feature of the single particle Hamiltonian (|5|) , which was completely ignored in Q , is that due to the 
weakness of disorder it almost conserves the momentum. It is natural to parameterize the single particle energy e by 
the momentum fc (again we are not interested in the very edges of the spectrum) 

e = -2icosfc . (9) 

Because of w <C t, the eigenfunction of (^|) ipe( n ) should be a linear combination of cos(fc e n) and sin(fc e n) with 
slowly varying amplitudes within the intervals small compared to h (JsJ) . Therefore, it is convenient to consider the 
plane wave basis 



Mn) = }2C(k)\k) , |fc) = — e *« , (10) 

where L is the total size of the sample (L ^S> h)- Now the amplitude C(k) has two narrow peaks (Afc ~ l/h) around 
k = ±k e . 

Being the Fourier transform of the localized oscillating function tp e (n), the amplitude C{k) should manifest some 
simple features. Let the function ip e (n) be localized around n = rig . Then one has 

C(-k) = C(k)* , 

C{k) = e tkn °C{k) , (11) 
C(k) = f((k-k e )h) , 

where f{x) is some smooth function concentrated around the region \x\ ~ 1. 

The ensemble averaged value of |C(fc)| 2 may be also extracted from the textbook |H| 



1 l\ dk 



1 de x " 



(12) 



where h is the Anderson localization length (||). This is the averaged value of |C(fc)| 2 and therefore it can not be used 
directly for our following calculations. However, (|l^) gives the proper estimate of the amplitude of the wave function 
in the momentum representation C(k) in the whole range of variation of the momentum k. 

In order to illustrate the physical origin of ( |l2| ) it is useful to rewrite the single particle Schrodinger equation in the 
plane wave basis: 

C W = ^«C( 9 ) , (13) 

<? 

where by w we have denoted the random potential (JsJ) and e is the single particle energy corresponding to this 
eigenfunction. In general this equation is not easier to solve than the original Schrodinger equation with the Anderson 
Hamiltonian (||). However, in order to extract the information we need, it is enough to observe, that the r.h.s. of ( |T^ ) 
is saturated by the only C(q) with q — ko ~ l/h, and (for negative q) q + ko ~ l/h > where fco is the first positive 
solution of the equation eo = —2t cos fco • Therefore, far away from the region \k ± fco| ~ l/h the amplitude of C(fc) 
is determined by the energy denominator in (13) in accordance with (|l2|). 



Now it is easy to write down the estimate for the amplitudes C(k). 



C(k)~\ Y if k±k ~l/h , (14) 
V L 



C(*)~^_L_ if l// 1 «|fc 0± fc|«l . (15) 

Everywhere except the Section 5 we suppose that L 3> h- The formula ([IJ) follows directly from the normalization 
condition ^ |C(fc)| 2 = 1. The denominator fc ± fc in ( p^| ) stands for e — e k in ( |l3| ) and all the other factors may be 
found by comparison of ( |l5|) with jl^ ) at |fc ± fc | ~ l/h- In general, for arbitrary |fc ± fc | ~ 1 



Now let us take into account the interaction Vmt ■ It is natural to use the basis of products of single particle 
localized states 



|l,2) = ^i(n 1 )^(n 2 ) . (17) 

This wave function describes distinguishable particles. The generalization to the case of bosons or electrons with 
opposite spins is straightforward. Moreover, in the following section we will effectively turn to the bosonic case 
because the Hubbard short-range interaction (Q) separates the symmetric and antisymmetric eigenstates and the 
Hamiltonian (Q) is invariant under particle permutation (m — > n-i , n-i — > n\). 

In terms of the summation over the original lattice-sites the matrix element between the states (0) takes the form 

(3,4|Knt|l,2)=l/5^^ 1 (n)^a(n)^3(n)^(n) . (18) 

n 

Each wave function here decays exponentially outside the segment of the lattice of length ~ l\. Therefore the matrix 
element (^) vanishes if not all the functions ipi(n) overlap. Everywhere in the following we consider only the matrix 
elements for overlapping states. 

In order to estimate the matrix element (p^), the author of Q assumed that each wave function "4>i{n) within its 
segment ~ l\ is completely random. Under this assumption one immediately finds the following estimate 

(3,4|^ nt |l,2) chaotic ~-3 7J . (19) 

On the other hand, the assumption itself about chaoticity of the wave function is evidently inconsistent with the 
accurate estimate of the Fourier transform (|l^) - ( |l6| ) . 

The difficulties in estimation of the matrix element ([l^) originate from the almost regular and fast oscillations of 
the single particle wave functions. It may be seen immediately from any toy example, that the matrix element for 
oscillating functions is usually much more suppressed than for random ones. On the other hand, for oscillating wave 
functions it is natural to consider the matrix element ( |l8| ) in terms of plane wave amplitudes ( [l0| ) 

(3,4|y i „ t |l,2) = ^^C 3 *(fc 3 )C^(fc 4 )C 1 (fc 1 )C 2 (fc 2 )^ 1+fe2 , fe3+fe4 . (20) 

ki 

Here the ^-function accounts for the conservation of the total momentum by the interaction Vi n t- One may use 
( |i"T| ) in order to demonstrate in terms of C(fc) that only the matrix elements between spatially overlapping states 
ipi, i>2, i>3, i>4 survive. 

Equation (|2(]) together with the plane wave amplitudes (|l4|)-(|l6|) is enough for the rough estimate of the matrix 
element of Vint- As we have said, we will consider only the overlapping states (no —no ~ h (|ll|)), otherwise the matrix 
elements decay exponentially. We have considered three different estimates (|l4|), (jl5|), and ( |l6| ) of the plane wave 
amplitude C(fc) for different ranges of variation of the momentum fc. Correspondingly, now we are going to present 
the three different estimates of the matrix element (pp|). First of all, consider the single particle states |1), |2), |3), |4) 
for which the momentum is almost conserved in the transition (|20|). This means that k 0l + k„ 2 w fc 03 + k 0i (the 
momentum k is connected with the single particle energy via (^|)). As we know each function Ci{k) consists of two 
narrow (Afc ~ l^ 1 ) peaks. Just for fc 0l + k 02 w fc 03 + k 0i all this peaks overlap and therefore 



"H k 02 k (>3 k 0l ~ l/l 
U fh\ 2 fL\ 3 U 



(341^112)--^^ ~ - . (21) 

Here the first factor U/L comes directly from (p0|), (li/L) 2 is the fourth power of the wave function at the maximum 
(|l4|), and (L/li) 3 is the effective number of terms in the sum in (po[). Suppose now that the momentum conservation 
is completely violated Afc ~ 1. In this case the peaks for all C(fc) could not overlap simultaneously. One of the 
functions should be taken at the tail (|ilf): 



ko 1 + k 02 k . 3 k 0i 



(341^112) ~ - ( i ) -(-] ~ . (22) 



5 



Here compared to (|2l]) we have replaced one of the plane wave amplitudes ( |l4| ) by (|l6|). The factor 4 in (22) should 
not be considered very seriously, it simply symbolizes that this replacement may be done in 4 ways. 

Finally, the most interesting case is if the momentum difference Afc = k 0l + fc , — k . 3 — k 04 is much larger than \jl\ 
but still is small compared to 1. To be more accurate one should define a few Afc-s , Afc = ±fc 0l ± k 02 ± k 03 ± k 0i . The 
matrix element will be enhanced as in ( p3| ) below if at least one Ak is much smaller than 1. We leave the consideration 
of this complication until Section 4. The formula (|l5| ) for the "short range tail" now allows one to find the following 
estimate 

l/h « |fc 0l + k 02 - k 03 - fc J < 1 , 
(3,41^,11,2) ~ ^ — — . (23) 

Just this last matrix element will lead to the enhancement of the two-particle localization length. We will return to 
the more accurate and detailed estimate of the matrix element for the case ( ^3| ) in Section 4. 

Note that all these estimates (|l|),@),(||) differ from those of § (3,4|F int |l, 2) - U/l\ /2 ©. 

The natural tool for the investigating the coherent propagation proposed in Q is the Thouless block picture [ p(i[ . 
Consider the two one-particle Id Hamiltonians with the interaction (|j) as the 2d-Hamiltonian on a large Lx L square. 
Because the Hubbard interaction (Q) affects only the diagonal n,\ — n% it is natural to use (for a moment) the center 
of mass variables \{n\ + TI2) and (ni — 712). Now the large 2d system should be mentally divided into square blocks 
of size lb x l), with lb ~ l±. Due to the exponential decay of the localized eigenstates at |n — n \ > l\ and the diagonal 
form of the interaction (0), it is enough to consider only one row of blocks with |m — 7T-2 1 < h/2. Thus our 2<i-system 
reduces to quasi-l<i one. The m-th block is defined by the inequalities 

|ni-n 2 |<| , (24) 

mlb < < (m + l)l b . 

About l\ r~j l\ states (^) fall down into each block. Thus the typical level separation within one block is A ~ t/l\. 
We suppose that the interaction (Q) mixes the states from the same block and with about the same amplitude the 
states from the nearest blocks. If one believes in the estimate (|l9|) Q| the derealization (enhancement of the coherent 
propagation length) follows immediately from the inequality 




(3, 4|Vi„ t |l, 2) c/iaot j C ~ Ayii > A . (25) 



However, unfortunately the matrix element ([19]) is inconsistent with the correct estimate. 

Consider now, what one can conclude from the correct estimates (^l"l)-(^3|)? In order to understand whether the 
enhancement exists or not, it is enough to fix one state ( p"^ ) belonging, say to the m-th block and to estimate the 
number of other states mixed to that one with amplitude ~ 1. First of all, consider the largest matrix element (pi]). 
This matrix element (34|Vi n t|12) ~ 1/Zi is l\ times larger than the typical level spacing within one block. Nevertheless, 
the momentum conservation necessary for ( f2l| ) reduces in the same l\ times the number of states in one block available 
for this transition. Thus 

A e//o ~ (3,41^*11,2) ~ i for |fc 0l + k„ 2 - k (>3 - k J < i , (26) 

h ' n 

which means that about one state is strongly mixed to the given one by this matrix element. 

So we are looking for the admixture to the state |1, 2) belonging to the m-th block. Let us divide all the ~ If states 
|3, 4) from the m-th and two nearest blocks into smaller portions in accordance with the momentum non-conservation 
Afc = fc 0l + k„ 2 — k„ 3 — k 0l . Into the n-th portion we will put the states with 2 n /l\ < |Afc| < 2 n+1 /l\. Due to ( p3] ) 
one may easily compare the effective level splitting A e f f n for each portion with the corresponding matrix element 

A e//l ~(3,4|V int |l,2)~-_ for i < |fc 0l + k 02 - k 0s - fc„ 4 1 < i , 

A e/A ~ (3,41^11,2) ~~ for £ < |fc 0l + k 02 - K 3 - fc 04 1 < 1 , (27) 



4/i h 



1 1 2 n— ^ 2 n 

A e //„ ~ (3,4|y int |l,2) ~ — - for — < |fe 0l + k f<2 - k„ 3 - fe J < — 



G 



The number of states falling into the n-th portion increases like 2™ but simultaneously the matrix element (|23|) 

decreases by the same factor. The number of rows in the formula ( |27|) is evidently ln(Zi) and therefore we see that 

each simple state |1, 2) ( |l7| ) may be mixed with ~ ln(ii) 3> 1 others. 

We will see in the Sections 5,6 how this ~ ln(Zi) enhancement may lead to the coherent propagation length (||,[|) 
The approach of ref. p2] is sometimes considered as one which allows to get the better understanding of the TIP- 

problem. Therefore before going further we would like to consider the validity of the method of (l2) in the light of 

our estimate for the two-particle matrix element (pl|)-(p^). 

The authors ref. [Q calculated the averaged Breit-Wigner width T for TIP in small L <C l\ sample. Connection 

with the full TIP problem is made via the relation 

h/h ~Tp , (28) 

where T is the Breit-Wigner width for the sample of the size L ~ l\ and p is the total density of two- particle states 
in the same sample. The authors do not derive the eq. (p8|), but find the support for it in the refs. jljME] [§)• The 
eq. (^8|) seems to be useful if the interaction would be able to mix with about the same probability all overlapping 
two-particle states. However, as we have shown in ( pT| , p2|j2^ ) the actual matrix elements between different states are 
of the very different magnitude. For example the total width Y in (^8|) is determined by a very small fraction ~ 
of all two-particle states. In its turn the main part of the density of states p is given by the states which may be 
mixed only by a very small matrix element ( f22| ) (in 1/y/h times smaller than is needed for the estimate of ref. jy). 

Thus, there seems to be no reason to join in one expression ( |2lj| ) the T and p which have so different physical origin. 
One may divide all the two-particle states into a classes so that the matrix elements from the given state 1 12) to 



any other within one class will be of the same order of magnitude. Just as we did in the eq. (27). Now the natural 
generalization of ( p8|) will be to introduce the large parameter (see (p7j)) 

g~ Yl FiPi-U/tMh) , (29) 

classes 

where Ti is the partial width for the state 1 12} to decay into the i-th class and pi is the density of states in this class. 
This g will work as the effective expansion parameter in the perturbative treatment of the TIP problem. It is however 
still a long way from eq. ( f^t ) to the accurate estimate of the coherent propagation length l^- This will be just the 
main aim of our Sections 5,6 to show that in order to find the TIP coherent propagation length one should most 
naturally exponentiate g and that h/h ~ exp(const x g) . 



III. THE STRONG COUPLING LIMIT 



The existing attempts to reduce the two particle problem to a random matrix one were all based on the assumption 
that the mixing of the simple states (17) due to the particle-particle interaction is sufficiently random. 

On the other hand, if there is no disorder (w n = 0) the exact eigenfunctions of the Hamiltonian (^) are easy to 
found in terms of the variables n\ — ri2 and n\ + ri2 . Among these solutions the eigenfunctions decaying exponentially 
at large \ni — ni\ form the molecular bound states sub-band. All remaining wave functions from the continuous 
spectrum will be formed by rather regular combinations of sin(ni j 2) and cos(ni.2), though having a finite kink at 
n i = n 2 ■ Thus for w n = the particle-particle interaction leads to considerable but rather regular rearrangement 
of the noninteracting product basis set (p"7|). It is evident that for finite but small disorder {\w n \ <C t) the regular 
modification of the wave-functions should survive, at least in some sense. 

This regular rearrangement of the noninteracting basis may be thought as the mixing of each state ( |l7j ) with many 
others via the largest matrix elements (|2|). The energy denominator for this many states in general is not small (up 
to AE ~ t) and thus each individual admixture is added to ( |l7j ) with a rather small amplitude. However, due to the 
large number of these effectively regular contributions, the final change of the wave function may be of the order of 
one. 

Thus we divide (although slightly arbitrarily) all the admixtures to the given state |1,2) into two classes. First, 
those with small energy denominator AE < These contribution are completely random as we show in the 

following section and lead to the large coherent propagation length. 

The second type of admixtures are those with large energy denominator AE > l//i. We see no mechanism, how 
these corrections may lead to the enhancement of the propagation length (although also we have no rigorous proof 
that such mechanism does not exist). Moreover, we are going to show in this section that these regular effects in the 
strong coupling limit should lead to the suppression of the two particle propagation. Up to now the authors of all 
the papers devoted to the TIP problem have examined only the monotonous dependence of I2 on U (say I2 ~ U 2 in 
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[|J and h ~ |C/| in ||). Among other papers the numerical result of ref. H may be considered as the most serious 
objection to our prediction ([3]). However, our main statement in this section is that all quantities characterizing the 
TIP problem should depend in some complicated way on U/t even for U ~ t. For example 7 in (0) for U < t may be 
of the form (we do not worry now about the sign of U) 



U fu\ 2 /U^ 3 



l{U) = li 1 +l2\jj +73 [jj +-.. , (30) 

where all 71,72, ... are of the order of one 7$ ~ 1. The same holds for the overall normalization of (Q). Therefore, 
because the authors of || have not taken into account this complicated dependence of I2 on U, we believe that they 
have used the irrelevant fitting function for the l\ dependence as well (see also the discussion [fL9||). 

Let us consider now the case of very strong Hubbard interaction (j^) U t. Consider also the new basis set of 
states relevant to this limiting case. First of all, the trivial subset of diagonal states evidently decouples for the large 
U 

\n) = §n,mSn,n 2 j (31) 

where n\, n,2 are the lattice variables for two particles. These states form a narrow molecular zone with energy close to 
U and the effective molecule mass close to ~ U/t 2 . If t 2 /U <§C w these particles do not move. If t 2 /U » w, Anderson 
localization takes place, though the localization length is small compared to 

The most surprising fact is that in the model under consideration all the remaining eigenstates in the strong coupling 
limit may be found exactly. Let ipi(n) be a set of localized single particle states. Then it is easy to construct the set 
of antisymmetric states 

|l,2) A = -^[Vi(n 1 )V2(n 2 )-^ 2 (ni)^i(n 2 )] , (32) 

which are obviously the exact eigenfunctions of the Hamiltonian (|j) even for arbitrary U. In addition to (|3l]) and 
( |32| ) let us consider the symmetric functions 

|l,2) s = -j=sign(ni - n 2 ) fyi ("1)^2(712) - ^{nxjipi (n 2 )] , (33) 

where sign(x) = x/\x\ and sign(0) = 0. First of all, the states (^l|)-(|33|) form the complete orthogonal set and are 
all localized. Any two vectors from (|l]) and (|3^) are orthogonal by construction. The states (|33|) are evidently 
orthogonal to ([52]) and (|l]). Finally, the integral of the product of any two symmetric states (|33[) coincides with 
the one for corresponding antisymmetric states. Therefore. (|3l|)-(|33]) may be used as a new basis instead of Jl7|)- 
Moreover, as we will show now the symmetric functions ( [331 ) are also the exact eigenfunctions of (Q) for U — > 00 (or 

u->t). 

The physical explanation, why (|3^) are the exact eigenfunctions for large U is almost trivial. Consider our two 
particles in Id as one particle on a 2d square L x L. For infinite U this particle simply could not penetrate through 
the barrier (Q) along the diagonal of the square. As we said before, the antisymmetric states ( |3^ ) are by construction 
the exact solutions. But now, because the two parts of the volume are separated by an opaque barrier, the symmetric 
states (||^), as well as any linear combination of ( |3^ ) and ([33]), become exact eigenfunctions of the Hamiltonian (Q). 

However, this simple explanation does not allow one to find the corrections to (^3|) for large but finite U . Therefore 
below we calculate the matrix element of the effective interaction between the two states |1, 2) s and |3, A) s for U t. 
Consider the total Hamiltonian H to t (f§) as a matrix in the basis (|3l])-(33). It is enough to consider only the symmetric 



states ( |3~l| ) and ([33|). The only surviving non-diagonal elements of this matrix are 

(n|# tot |l, 2) s = -V2t[tyi(n + 1) - - l))^(n) - (V> 2 (n + 1) - ~ l))^i(»)] • (34) 

Besides that there are evidently only the diagonal elements 

(n\H tot \n) = U + w„ » U , (35) 
(3,4| s i7 tot |l,2) s = ei +e 2 . 

Thus the mixing of the two states |l,2) s and |3,4) s may appear only in the second order of perturbation theory. 
Consider for simplicity the case e x + e 2 ~ e 3 + e 4 . The simple perturbative calculation immediately leads to the 
effective interaction 
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(3,4| U eff \l,2) =^ [r - — = _____ (36) 



Ci(k 1 )C2(k 2 )C^(k 3 )Cl(k i )[sink 3 — smk 4 ][smk 1 — sinfc 2 ] 



This effective matrix element should be compared wit h (P0|) . First of all we see that the strength of the interaction 
U has been replaced by t 2 /U- Thus the corrections to (|33[) are of the order of t/U -C 1. In particular the estimate of 
the matrix element (|2l| - p3|) evidently holds for ([36]) up to the trivial replacement U — > t 2 /U (see also ((H) below). 

Another important difference between ( pll| ) and ( |36| ) are the sin(fcj)-s, which can not lead to any enhancement and 
even tends to suppress the matrix element for e l ~ e 2 and e 3 w e 4 . 

Thus we see that the problem of TIP should be treated within the two sufficiently different (somehow dual) 
approaches in the weak ( [l8|]2C| ) and strong ( |36|) interaction cases. However the l\ dependence of the matrix element 
estimated in ( f2lf|23| ) coincides for both limits. 

In the intermediate region U ~ t both formulas (^0|) and ( |36"1 ) should be modified. In order to illustrate the nature of 
this modification consider for example the case U < t. Let us consider also the eigenstates with the two-particle energy 
close to some value E$. It is convenient to consider the narrow strip in energy of the width 2Se ~ t/l\ around Eq. As 
we have told at the beginning of this section the states within the strip \Eij —Eq\ < Se will be mixed chaotically. The 
only role played by the other states \Eij — Eq\ > Se is to renormalize regularly the effective matrix element between 
chaotically mixed ones. Let the states |1, 2) and |3, 4) belong to the strip \E%2 — Eq\, I-E34 — Eq\ < Se 2) and |3, 4) 
are the simple products of noninteracting localized single particle states). Than one may write perturbatively the 
renormalized interaction 

(3,4|V ren |l,2) = (3,4|V int |l,2) + (37) 
+ y (3,4\Vi n t\i,j)(i,j\Vi n t\l,2) 

\e^e \>s e E °- E « 

The most interesting for us situation is when the momentum in the transition |1,2) — > |3,4) is not conserved and 
the leading contribution in ( |37| ) is described by ( p3| ) . In this case one of the matrix elements in the sum in (^) may 



be large like U/h (21) (and the small factor in it may be compensated due to a number of terms in the sum) 
However, as one may easily see, the second matrix element in the sum should again be of the form (p3|). Thus we see 
no reason, why the renormalization of the bare matrix element (|l^,^,^) will not be described by the series of the 
kind 1 + (U/t) + (U/t) 2 + ... . 

It is clear that the states ([33]) should be also the proper solutions for the continuous version of the model ([|) with 
strong repulsive interaction Vint — G5(x — y). Of course in this case the energy e x +e 2 should be also sufficiently small 
(e, + e 2 <C G 2 ). In fact the continuous model may be simply considered as the low energy limit of our lattice model 



To conclude this section let us say again that the comparison of our eq. (|36|) with ([20|) shows that the crossover 
from weak to strong coupling regime takes place somewhere at U/t ~ 1. The authors of refs. |i 12 have mentioned 



that they would expect some deviations from their predictions at large U. However as we can understand, they have 
in mind just the case \U\ ~^>t. Moreover, one would even hardly found from |||l2|], what is the sign of these expected 
corrections. 



IV. ESTIMATION OF THE MATRIX ELEMENT 

The rough estimate of the matrix element ( p3| ) in principle is enough to formulate the effective matrix model which 
will be considered in the following two sections. Nevertheless, because the distribution of the matrix elements is of 
crucial importance for our result, wc would like to present in this section a more accurate estimate than (p3|). We will 
consider both the matrix element of the interaction potential (Q) in the basis of the simple products of noninteracting 
single-particle states ( p7| ) and the matrix element of the effective interaction in the strong coupling case (|36|). 

Although, our formulas ([53]) , (|54"|) below are too detailed to be used for the mapping onto the matrix model (Sections 
5,6), they may be a good base for the comparison with the numerical experiment, say, of the kind mentioned in 0. 

Our main problem is to make a reliable estimate of the matrix element with oscillating wave functions. As we have 
said in Section 2, we are mostly interested in the case when the four plane wave amplitudes Ci{k) do not overlap 
strongly in the matrix element ( p0| ) . This means that at least one of the wave functions should be taken at the tail 
(|l5|). For this case the Schrodinger equation ([l3]) may be further simplified 







C(k) tail = J2 



(k\w\q) 



C(q) c 



(k\w\q) 



\q—k \<S 



|g+fco|<<* 



C{q)c 



(38) 



where \k ± k \ 3> 1/h and S ~ l//i (say <5 = 5/Zi ). Here we have divided each wave function C(fc) into two parts. 
One part is "the core", or the largest part of the wave function coming from \k — k \ < S and |fe + k \ < S. Another 
part is the tail at \k ± A; | > 6 . 

The matrix element of the random potential between the two plane waves which appeared in (pq) is equal to 



(k'\w\k) = (0\w\k-k') 



L 



(39) 



For the quantitative calculation of the matrix element of the interaction potential it is important that (p9[ ) depends 
only on the one variable k — k! . The distribution of (k'\w\k) is an usual normal one with variance (see (H)) 



\(k'\w\k)\ 2 = w 2 /L. 



(40) 



es : 



Formula (|38|) allows us to demonstrate one important difference between the "core" and the "tail" of single particle 
wave function. The main part of the wave function is almost regular (one should not worry about the trivial randomness 
due to the position of the center of the wave packet no ( |ll| ) because we consider only overlapping functions). On the 
other hand the tail ( |3^ ) due to (0|u)|fc — k') is sufficiently random. More accurately, C(k) and C(fc') are correlated 
only for k — k' ~ X/l\ and are not correlated for \k — k'\ 3> 1/h- 

This observation allows us to define a relatively simple procedure of averaging the squared matrix element 
The mathematically rigorous way of averaging the |(3,4|Vi„t|l, 2)| 2 would be to introduce eight collective variabl* 
centers of the wave packets n 0l , n 02 , n 03 , n 0i and the energies e 1 ,e 2 ,e 3 ,e 4 . The "rigorous" averaging should then be 
performed over the ensemble of random Hamiltonians (||) for fixed values of these collective variables. However, the 
center of the packet and the energy are almost completely determined by the smooth core of the wave function. The 
actual randomness of the matrix element, as we will see below, is hidden in the weak tail ( |38| ) which almost does not 
affect no and e. Therefore in this section we are going to consider the distribution of the matrix elements (3, 4|Vi„t|l, 2) 
for "frozen" cores. In terms of the single particle matrix elements ([59]), it means that we have frozen the values of 
(0|u)|A) in a few small (~ 1/Zi) segments around A = and A = ±2fc . (i = 1 — 4) and have performed the averaging 
over realizations of all the others elements (0|w|A'). 

As in the previous sections, consider the transition between the states |1,2) and |3,4). Let now |1,2) and |3,4) be 
the symmetrized simple products of the kind ( |l7| ) : 



1 



1 1,2) = -^=(^1(^1)^2(^2) + ^2(^1)^1(^2)) 



(41) 



This symmetrization simply adds an overall factor 2 to the matrix element (18 2C|). 

We are interested in the matrix elements between the approximate eigenstatcs which may be mixed strongly by the 
interaction. Therefore, because the largest matrix element ( pl| ) is of the order of ~ l/h, the product states should 
be almost degenerate £! + e 2 — e 3 — e 4 ~ tjl\. Moreover, because we have seen in Section 2 that the smallest matrix 
elements ( p2| ) themselves could not lead to any considerable enhancement, it is enough to choose |e x — e 3 | -C t and 
|e 2 - e 4 | < t (although still t/h < h - e 3 | Je 2 - e 4 |). 

We have divided each wave function in ( p^ ) into two parts: a strong "core" and a weak (~ w/Ae) "tail". Due to 
momentum conservation in the matrix element (pfj[), direct transition between cores is not allowed now. Equation 
( |38| ) shows however that this transition effectively may take place in the first order of perturbation theory in (k\w\q). 
Substitution of (J3q) into (Effl) (with the additional factor of 2 due to the symmetrization (fidh) gives 



U 



(3,4|F mt |l,2) =2-£ CabMiqJCiiqJCMS^ 



3+94-91-92 



(42) 



x(u w 



|A) 



A qi—core 

1 



t 9l +A 



Here the summation over qi includes only the largest plane wave amplitudes (|Tj) C(q) ~ \JT\JL. Each of the four 
terms in the figure brackets in ((42|) corresponds to the replacement of one of the C(q) in (EOT) by the perturbative 
formula (38). However, due to (|39|) all these four contributions have the same random weight(0|w| A). 

In Section 2 we have classified all the matrix elements according to the value of momentum non-conservation in 
the transition Ak. In order to make the simplest estimates up to now we made almost no difference between the 
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momentum k and the single particle energy e = — 2tcos(k). However, the matrix element ( |23| ) turns out to be 
enhanced just due to the small energy denominator. Each wave function in momentum representation consists of two 
peaks k ps ko and k ps — k a (we choose k to be positive). Thus, effectively in the equation in addition to the 
four terms in the figure brackets, we have 4 x 4 = 16 transitions 



i k 02 \ 
~ko 13 kv 2 I 

"'Oi 1 "-02 I 



/ ^0 3 i ^0 4 



V k 0l , fc. 



"2 



V-fc, 



^03 ! "-04 
3 7 ~^0 4 



(43) 



between pairs of narrow peaks (cores). The transitions between these pairs became possible due to perturbation theory 
(|38|), but not all of them (although still many) will be enhanced by the small energy denominators. On the other 
hand, the perturbative formula (|38| ) also contains the random matrix element (0|u)|A). This A will be different for 
different transitions from ( f43| ) and thus different contributions will not interfere in the averaged value of the squared 
matrix clement of Vi n t- This means that one may calculate separately the contribution of each two- by- two pair in 
(||) into |(3,4|V/ mt |l,2)| 2 . 

Consider first the case (k 0l , k 02 — > k 03 , fc 04 ) (only this case was supposed to be considered in the rough estimate (|2 



of Section 2). All the energy denominators in 
this means that 



now are small compared to the typical energy ~ t. For example, 



-ei+A 



(fc 0l - <?i - A)w 1 w -A«! ps -(&„ 3 + k (>4 - k 0l - fc Jui 



where we have introduced the velocity (see 



Substitution of i 



de 
dk 



= 2t sin k 



(44) 



(45) 



|) into (|42|) for the transition (k 0l , k . 

U ^ 



(3,4|VW|1,2) ++ , 



E 



k o 3 ,k J gives 
C, 3(93)C , 4(g4)Ci(q' 1 )C'2(q 2 )5A,«, 3 +g 4 - gi - q2 



(0|«>|A) 



kn 1 



^0 2 



ko 3 k 04 ^ t>i 



1 



(46) 



Here (k 0l + k 02 — k 03 — k 0i )~ 1 is the large factor considered in eq. ( ^3*1) of Section 2. The subscript ++,++ stands for 
the signs of (fc 0l , k (>2 ) and (fc 03 , k 0i ) in the transition. 

However, if k 0l + k 02 — k 03 — k 0i <C 1 and e 4 — e 3 pa e 4 — e 2 -C i one has 



Vi - v 3 



v o 3 



ko A 



(47) 



Thus the small denominator in ([16|) turns out to be compensated by the strong cancelation of the terms in the 
figure brackets. On the other hand, this cancelation is specific only to the transitions (k 0ll k 02 — > k 0;i ,k„ 4 ) and 



( k Ql , k 02 



v o 3 > 



-fe D4 ) (or (++,++) and 



-)). For example, let us consider the case (k 0l ,k 02 ) — > (fc 03 ,— fc„ 4 ). 



Now A in ([42[) is of the form 



A 4 



— ^o 3 "»o 4 ^Oj 



"2 



(48) 



Here A ++ + _ is not small: A ~ 1. However, for A ++ + _ two of the energy denominators in ( fl2|) are of the order of t 
but other two are still small: 



e<j 2 +A 
e 94 -A 



e — fco 2 e feo 2 +A++, + - 

e *-'o 4 — e -*;o 4 -A ++i+ _ 



(fco 2 - fc o 4 + fco 3 + fc 0l )«2 
(M) 4 ~ fco 2 + ^o 3 — k 0l )l!4 



(49) 



Simple substitution of ( |49| ) into (42) gives 

(3,4|^ mt |l,2) ++ , + _ 



'xE E E 

A <?l,2,3~fco li2 , 3 94~-fco 4 



,93+94-gi -g 2 
2(k 0l - k 03 )(0\w\A) 
v 2 [(k 0l — k ( , 3 ) 2 — (fc„ 2 — fc r , 4 ) 2 ] 



(50) 
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where we have made use of v 2 ~ v 4 . Now it is convenient to go back from momentum to coordinate space. Equations 
(|Io|) , (|39|) allow one to rewrite (|5(]) in the form 

(3,41^*11,2) (51) 
iU v t v 2 



,w n ^i(n)ip2(n)ip 3 (n)ipj (n) 



Here we have used: e x — e 3 » e 2 — e 4 and, e.g., ei — e 3 s» t>i(fc 0l — ^o 3 )- V ;+ an d V' - arc the positive and negative 
frequency part of the single particle wave function respectively (ip + = ip~*)- 

As we have mentioned before, due to the random factor (0|u)|A), the contribution (pQ) docs not interfere with 
others in the squared matrix element. On the other hand, due to the random factor w n in (|5l|), there are no problems 
now with cancelation of the sum of regular oscillating terms. Thus, V' =t m ( pl| ) are the almost regular solutions of the 
single particle Schrodinger equation ( "cores" ) which as we have seen are responsible mainly for the global features of 
the eigen- vectors such as, e.g., the energy and the position of the wave packet. Therefore one should naturally average 
the squared matrix element via 



lX>„V^(n)^(nWr(n)Vtf(n)| 2 = (52) 

n 

l^^(n) 2 ^(n) 2 V3(n) 2 V4(n) 2 . 



w 2 

16 



We have considered in (|5l|),(p2|) only the contribution of the transition with signature (++,H — ). It is easy to 
show that identical contributions to the variance of the matrix element come from the transitions ( — , — h), (h — ,++), 
(h — ,--), (++, — h), ( — ,H — ), ( — h++), and (h — ,--). Finally, all these eight contributions give 



|(3,4|y int |l,2)P = - ? 2 1 \ \ MnYMnrMnYMnr (53) 

As we have said before the averaging here (the intermediate averaging) is performed over the Fourier components of 
the random potential (0|u>|A), which do not contribute sufficiently to the "core" part of the wave function, namely 
ii| A| ^> 1 and l\ |A ± 2k . \ ^> 1. The equation ( |53] ) may be used if t/l\ <C e 1 — e 3 <C t but the energy non-conservation 
in the transition is very small e t + e 2 — e 3 — e 4 ~ tjl\. In particular one should not use this equation for calculation of 
the Breit-Wigner width T considered in jl2| . As it was said at the end of Section 2 the total width T is determined 
by only the few largest matrix elements. In our notations these are the matrix elements with e 1 — e 3 ~ t/l\ and (|5^ ) 
may be considered in this case at best as only the order of magnitude estimate. 

One may easily repeat the calculation (50)-(j53|) for the less restrictive case e x — e 3 ~ e 4 — e 2 ~ e x + e 2 — e 3 — e 4 <C t. 



|(3,4|^„ t |l,2)|2= (54) 

4t/W - - .. (Cl = ± (£2 ' U) \ 2 - , J2MnYMnYMnYMnY 
v 2 v 2 ((e 1 - e 3 ) 2 /v 2 - (e 2 - e 4 ) 2 /« 2 2 ) n 

In the following two sections we are going to consider the mappings of the two-particle problem onto some special 
random matrix models. However, these matrix models, although are expected to reproduce qualitatively the main 
features of the two particle cigenstates, certainly could not be used for any quantitative calculation. Therefore it seems 
rather probable that one still should perform in the future the accurate calculation of I2 with the original Hamiltonian 
(Q) . The numerical verification of our expression (|||) seems to be a good starting point for such calculation. 

In a similar way one may find the variance of the matrix elements of the effective interaction ( p6| ) in the strong 
coupling basis (|3^) . We skip the details of this calculation and give only the result 

At 4 w 2 



|(3,4|^ mt |l,2)sp = ^ 2(ei _ £3)2 \: ^MnrMnrMnYMnr ■ (55) 

Here we have neglected e 1 + e 2 , compared to U in the denominator. 

The matrix element in the simple basis of noninteracting states (p0|) , (|53"|) is proportional to (1^ — u 2 ) _1 . At first 
glance this factor may effectively strengthen the interaction for e t « e a . However, in the strong coupling limit ( |55"| ) all 
denominators (Wi — v 2 )~ 1 disappear. The reason for this more regular behavior of the matrix element in the strong 
coupling limit is obvious. If v ± is close to v 2 , the relative kinetic energy of the two particles is small compared to the 
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interaction (fjf) which is effectively equivalent to the strong coupling case. In Section 3 we have introduced the new 
basis se t (plfl33| ), which for U ^> t is relevant for all two-particle states. On the other hand, the factor (v t — Va)" 1 
in (pp|) , (|53|) indicates that even for t ^> U but e, sa e 2 one should modify the product basis (Ell). This modification 
consists of the strong mixing of many states ( |ll|) with e± w e 2 in order to get something like (p3|) . However, we do 
not see how this mixing (rather regular in fact) may lead to the enhancement of the coherent propagation length. 

The author of || , in addition to the particle-particle propagation has proposed to consider the coherent propagation 
of the particle-hole pair. For a particle and a hole with excitation energy small compared to Ef one naturally has 
Vi ~ v 2 . However our strong coupling basis can not be directly applied to this case because of the wave functions 
(^T|)-(|33|) do not describe the simple particle- hole excitation above some vacuum. Thus the coherent propagation of 
particle and hole may be very different from the coherent propagation of two particles. 



V. EFFECTIVE MATRIX MODEL. SHORT SAMPLE 



The block picture considered at the end of Section 2 allows one to perform easily the mapping of the two particle 
Hamiltonian (^) onto the corresponding random matrix model. However, almost no results may be obtained analyti- 
cally for this matrix model as we will see in Section 6. Therefore, we still would like to consider a simplified version of 
the problem of coherent propagation. Namely, let two interacting particles move inside the sufficiently short sample 
L ~ lx- In this case all particles are free to move throughout the whole sample and there is no room for the interaction 
induced derealization. However, one may be interested in the chaotic mixing of many simple states (like ( p"7j ) or (|3^)) 
with different single particle energies e. The natural quantity, characterizing the complexity of such states, is the so 
called inverse participation ratio - U pr . The explicit definition of U pr reads: 



^ = EKI 4 . ( 56 ) 

where a n are the normalized coefficients of the expansion of the chaotic eigenstate over the noninteracting two-particle 
states. 

In order to find U pr we are going to replace the solution of the exact Schrodinger equation (|]) by the relevant 
random matrix problem. First of all, as we have seen before ( pl| ) , the matrix elements of Vmt never exceed the 
value (3, 4|Vi„t|l,2) ~ U/lx (more concretely: max(3, 4|V e //|l, 2) ~ min(t/, t 2 /U)/li). Therefore we may restrict our 
attention to the consideration of only the part of the simple states, having their total energy within the narrow band 

Eq - 8 E < e 1 + e 2 < E a + S E , (57) 

where Se ~ t/l\ . The total number of such states now is N ~ lx ■ The matrix elements of the interaction potential Jig ) 
or (|3^) between these states form the random matrix. The main result of our previous consideration (p3|) , (p3|) , (po| , 
is that the statistics of the elements of this matrix is not uniform. In order to demonstrate explicitly the structure of 
the interaction matrix, one should simply enumerate the states ( |57| ) by assigning each one a number n in accordance 
with, say, the energy of the first particle: n^) > n'(e[) if e t > e[. Now, due to (|23|) , (|53"|) , (|55|) , the non-diagonal 
elements of our matrix M will be Gaussian distributed with second moment 



i^j ; MijM kn = (SjkSin + Sj„Sik) ■ (58) 

\ l J) 

In general u ~ 1/lx ■ In the weak coupling limit (U <C t) one has u ~ Ul^ 1 ( |53| ) and in the strong coupling limit 
(U ^> t) one has u ~ t 2 !/^ 1 ^ 1 (|55|). We have omitted the slow dependence of u on i + j ( |5^ , |55| ). This complication, 
if taken into account , should not change our main conclusions. The diagonal elements of the matrix V are uniformly 
distributed within the range 

M nn = e„ , -6 E < £ n < $e , (59) 
Se = Gu . 

Here we have subtracted from the total energy of each state the trivial constant E (|57j). We have also introduced in 
(59) the dimensionless parameter G = Se/u which allows to compare the strength of the diagonal and non-diagonal 
(58) elements of the matrix M. Analytical results in this section (as well as the poor analytical estimates of the 
following one) will be found only for G > 1, Moreover, just for G ^> 1 the correspondence between G and the 
parameters of the original Hamiltonian (Q-^) is most transparent. Due to the duality between weak ( |20| ) and strong 
(|36|) coupling cases G turns out to be the double-valued function of the strength of interaction U 
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G = const— , for weak coupling U <C t 
G = const^- , for strong coupling t <C U 



(60) 



The numerical constant in both cases is of the order of 1. The crossover from weak to strong coupling regime takes 
place at some G ~ 1. Thus the matrix model (58 5^) with very small G does not correspond to any physical limit 
of the Hamiltonian (ji[^) . In the crossover region instead of the matrix elements (|20| ) and ( |36| ) one should use some 
renormalized new formula for the effective interaction like, e.g., (37). This effective interaction should account for the 
strong (~ 100% for U ~ t) renormalization of both eqs.(|20|) and (36) due to the regular mixing with many two particle 

Namely these are the two-particle 



1 corresponds to a phase 
fl), the actual small parameter 



states which were not included into consideration by the matrix model 
states with \e 1 + e 2 — Eq\ > Se ■ However this renormalized interaction will naturally lead to the same effective matrix 
model (|5^,|5^) with G ~ 1 (see also ( |37| ) and discussion around). Also G is a function of only U and t and at least for 
w -C U, t we do not expect any considerable dependence of G on the strength of disorder. 

Matrices of the kind (58) are called Power-law Random Band Matrices (PRBM). Random Band Matrices have 
been the subject of intensive investigation during the last years ( for review see p2|). Power-law matrices with 
Mij ~ — j\ a also may be considered as a special type of random band matrices. However our matrices (^8|) with 
Mij ~ — j\ play an outstanding role among other PRBM because just the value a 
transition from localized to delocalized eigenvectors |Q . 

As we have tried to show in Section 2, due to the hierarchy of the matrix elements (t 
in our problem, e.g., for the weak coupling case turns out to be not U/t but (U/t) ln(/ij. The aim of this section will 
be to find analytically k pr (|5^) for the case (U/t) ln(Zi) ^s> 1 ^s> U/t (or equivalently ln(Zi) ^> G ^> 1). The method 
used for this calculation follows the Renormalization Group approach of ref. fl^| . 

Thus we consider now the case u <C Se ~ t/h ■ At a first glance, in this case one has to treat the non-diagonal part 
of the matrix M (pq) by perturbation theory. However the simple perturbative estimate gives a very small correction 
hpr — 1 ~ (u/Se) ^to the inverse participation ratio. On the other hand, no matter how weak is the interaction 
u , with small ~ U/Se probability the interval between the two diagonal elements of M may happen to be of the 
same order of magnitude as the non-diagonal one: e n ~ £fc ~ M n k- These rare events, which are not described by 
the usual perturbation theory, turns out to be the most important for the chaotic mixing of original non-perturbed 
eigen- vectors. Because the density of strongly interacting pairs of levels is very low for u <C Se, one may neglect in 
the zeroth approximation the effect of triple interaction, as well as the more complicated events. Thus, it is enough to 
consider the 2x2 part M of the full matrix M with diagonal matrix elements distributed uniformly between — Se and 
+Se, and the Gaussian distributed non-diagonal elements Mi,2 = -M2,i = v with the mean squared variance a 2 
. The eigenvectors for this matrix are easy to find 



M 



E\ v 
v e 2 



M 



Oil 
«2 



A 



Oil 

a 2 



2v 2 



1 



(61) 



The second solution is obtained by the permutation a\ 
the eigenstate (^lj) is 



y/(e!-e 2 ) 2 +iv 2 sfisx - e 2 ) 2 + ±v 2 T (ei - e 2 ) 

2 , a 2 — > — ct'x ■ The participation ratio corresponding to 



a 



= 1 - 



2v 2 



(si 



£2 



Av 2 



(62) 



It is easy to find the averaged value: 



P 



dv 



2vrcr 



exp 



v 
'2a 2 
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7r a 
2fe 



(63) 



This effect turns out to be much more important than the naive perturbation theory estimate 1 — P ~ (a /Se) 2 ■ 

Let us return to the N x N matrices ( p8| , [59| ) . In (^3|) we have found the averaged participation ratio for each of 
the two eigenvectors of M ( |6l| ) caused by the non-diagonal element v. Now each of the non-diagonal elements of the 
full matrix M m k leads to the same correction to the participation ratios of two from the TV eigenvectors up to trivial 
replacement a = u/\m — k\. Therefore (G = S E /u @) 



Pa = 1 — 



7T 2 1 

2 G ^ « 



l-^ln(Af) , iln(AT)<l 



(64) 
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Because of N ~ h, we see again that our actual expansion parameter is U/t ln(Zi) (weak coupling) in accordance with 
the prediction of Section 2. 

Physically the calculation of the participation ratio for our matrix model resembles the calculation of the Partition 
function of the interacting rare gas. The subscript 2 in P 2 in (|64| ) means that we have taken into account only one 
double event e n — £k ~ M n k and have completely ignored all the triple collisions, four-particles collisions, and so on. 
That is why the inequality describing the range of validity of ( |64| ) is so restrictive. However, let us consider now the 
case when only G _1 is small, but not ln(7i)/G. Consider the triple event. Suppose that the three diagonal elements 
Ei,ej,£k are close enough, so that the basis vectors e^e.,, e& may be mixed strongly by the interaction. Let us also 
have for defmiteness \i — j\ < \i — k\, \j — k\. The crucial observation, which in fact allows to find li pr is that, due to 
G> 1, and the decrease of the second moments (|58|) in the large ln(ii) limit one may strengthen the inequality: 



\i-j\<\i-k\,\j-k\ -> \i-j\<S:\i-k\,\j-k\ 



(65) 



Now the calculation of the triple contribution to P may be divided into two stages. First, the vectors and ej are 
mixed by the rather strong non-diagonal element My and form the two new basis vectors e.\ and e 2 in accordance 
with (|l]) 



§1 = ot i e l + ajej , 
§2 = — oijei + otiej . 

After that, one of these states (say ei) is mixed with by the much smaller elements Mjfc and Mj&. Due to (|5 
( |65| ) , the variance of the effective matrix element describing this mixing reads 

• 1 I ~ ~ 



|(1|M|£;)| 2 = \oiMa + ajMj h \* = 



(i-k) 2 {j-k) 2 (i-k) 2 {j-kf 



(66) 



and 



(67) 



Thus we see that the secondary interaction in the two stage triple event ( p5| ) effectively is not affected by the first 
interaction. 

This picture with the decoupling of the secondary interactions shows the applicability of the Renormalization Group 
approach for the calculation of U pr . Our aim in fact is to diagonalize the matrix M. Let us divide this procedure into 
N steps. At the first step let us diagonalize the three diagonal matrix e n , M nyTl+ i, M n<n —\. Because of e n ^> M nyTl ±i 
in general this procedure of diagonalization reduces (approximately) to the simple nullification of the sub-diagonals 
M n n ±i. Only for small part of states the diagonalization is nontrivial. Suppose, that e, — Ej+j ~ M i i+ i for some i. 
Now for these two states i and i+1 one has to use the exact procedure (|6l]). This diagonalization of the 2x2 sub-matrix 
should be followed by the orthogonal transformation of the remaining part of the matrix M : M n ^>n+i 1 Mn,k<n-i- 
Namely one should perform the following transformation of the two rows and two columns 



M„ 



M, 
M, 



i-\-l./n 



M, 



n,i-\-l 



axMi, n + 
-a 2 Mi n 



-I- aiM i+ i >n 



(68) 



OL\M n i + a 2 M,, 



-a.%M n i + aiM, 



where the index n enumerates the elements of the z-th (i + 1-th) row or column. ot\ and a 2 here are the same as in 
( |6l| ) |ai| 2 + \a 2 \ 2 — 1. However, due to (|65[), the statistics of the matrix elements M n ^> n +i, M n ^< n -i is not changed 
under the transformation (see (|67[)). In a similar way one may diagonalize out the third diagonal M njn ±2, the fourth 
diagonal M n>n ±3, and so on. 

Let us find, how the averaged participation ratio P flows under such gradual diagonalization. It is convenient to 
write the participation ratios of all N states in one column {p}. Now the gradual diagonalization of the matrix is 
accompanied by the gradual change of {p}, which we illustrate by the formula: 
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(69) 
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where the p k 's are the participation ratios for the sufficiently rare "important" events happening at the fc-th step [g3| . 
The probability to have the "important" event of the amplitude p is described by ( |62| , |63|) and the p k at each step are 
uniformly distributed over the column {p}. Combining together (M) and (p9), one finds: 




p=-^ = f[(l-^)~ N'^pr ~ I, lipr , (70) 

Hpr 



Of course this result was found for G ^> 1. Nevertheless the number of chaotically mixed states due to (|70|) may be 

very large U pr 3> 1. Most probably for G not large the formula U pr ~ l x lvr would be still valid, although 7^ pr in this 
case will depend on G in a more complicated way than in (|70|) . At least our numerical results which will be considered 
in the following section support this conjecture. It is also possible to extend the renormalization group approach of 
Jl7j used in this section in order to look for 7,^ (|70f) as a series in 1/G. 



VI. EFFECTIVE MATRIX MODEL. LARGE SAMPLE 



Now at last we are going to consider our main problem: the coherent propagation of two interacting particles in 
a very large disordered sample. In order to get the adequate matrix model we have to combine the block picture 
considered at the end of Section 2 with our understanding of the hierarchy of the matrix elements (|23|) , (^3|) , (|5q) . Any 
pair of overlapping single-particle wave-functions may be associated with one block from the row of blocks described 
in (^J). As in the case of the short sample, considered in the previous section, it seems to be enough to take into 
account only the chaotic mixing of the states from the narrow zone (|57|). Thus again there are effectively ~ l\ states 
in each block. It is natural to take into account only the interaction of the pairs belonging to the neighboring blocks. 
Therefore the exact Hamiltonian (Q) may be replaced by the block-wise random matrix 



/ M%, Ml, , 
Ml, Ml Ml 



M 



\ 



, Ml, Ml, Ml,. 
0, 0,M 3 b ,M 4 b ,. 



/ 



(71) 



Here M\ and M\ are the random If, x 4 (note that lb ~ h) matrices describing the transitions within i-th and between 
i-th and i + 1-th block respectively. In general the total size of the matrix (~ L) is supposed to be much larger than 
all the other lengths which may appear in our problem. 

In order to exploit the hierarchy of the matrix elements of the interaction one has to enumerate properly the simple 
states within one block as it was done for the short sample in the previous section: n(ei) > n'(e\ ) if e± > e^. After 
that statistics of the diagonal block matrices M\ will be again described by the formulas (^8|),(^). Moreover, the 

non-diagonal square blocks M\ from ( |7l| ) will be also described by ( |58| ) (of course the non-diagonal blocks M have 
no strong diagonal (|59|)). 



It is convenient to "smooth out" the block- wise matrix (71) so that 



< < 2l b , M tJ M kn 

where F(i — j) = 



G <e n <G , 

(SjkSin + 5 jn 8 ik )F(i - j) , 
a 1 



(72) 



{i-jY (i-j-h) 2 (i-j + h) 2 



2/ 6 < , M K 







Here we use the dimensionless parameter G introduced in (J5^). The connection between G and the parameters of the 
original Hamiltonian (^,0) is described by (|60|). In order to go from ( |7l| ) to ((7^) we have changed only a small fraction 
of the matrix elements which in the large logarithmic case ln(/i) ^> 1 are not important for the coherent propagation 
length. The amplitude a of the central peak of F in ( ^2[ ) is of the order of one: a ~ 1 (a > 1). 

The elements of the matrix ( |72| ) vanish, if \i— j\ > const x 1%. Therefore, as for usual band matrices its eigen- vectors 
should have some finite localization length. Moreover, because of the size of the block in ([7l|) is li, ~ l\ this localization 



1G 



length measured in vector indices should be proportional to the coherent propagation length I2 measured in the sites 
of original lattice (^). Thus the quasi-ld matrix Hamiltonian ( f72[) accounts for both the dimension of our sample 
and the complicated structure of the matrix elements of interaction (pjj) , (|53|) , j55|) . The central peak of the function 
F(i — j) (|7!) has the same form as the second moment (^8|) for the matrix model considered before for the single-block 
sample. However, as we will see below just the satellite peaks ~ (i — j ± 4) 2 in (72) are mainly responsible for the 
large localization length. Therefore, for future discussion and for numerical simulations wc have further simplified the 
model (72), so that 



M nn — £ n 
i < j < i + l\ 
i + h<3 



-G < e n < G , 

~r — t~? &jk3in ~\~ &jn$ik 



(73) 



U 



h + l) 



Mi. 



, 



Mi. 



This matrix model looks quite similar with that considered in the previous section (p8|) , (|59|) . The main difference is 
that in (58) the elements of the matrix decrease in the power law fashion as one goes farther from the diagonal, while 
in ( |73|) the matrix elements reaches the maximal amplitude at the two sub-diagonals j = i ± l\ and decrease like the 
power of \i — j ± l±\ towards the main diagonal. The model ( |73| ) should reproduce all the important features of the 
( [72] ) and thus will lead to qualitatively correct description of coherent propagation. 

The band matrices with Cauchy distribution of the matrix elements were also considered in [ B4[. This model due to 
the power law tail of the probability to find the large matrix element looks much closer to our (|75|) than the usually 
considered SBRM matrix models However, as is confirmed also by the author of [Q, for a rigorous application 

to the TIP problem more rigorous investigations are required. 

For G»l one may try to use the approach of Section 5 in order to find at least the inverse participation ratio li pr 
for the model (f73f). Moreover, if ln(h) <C G only the rare double interactions of the energy levels should be taken into 
account and the models (|5^ ) and ((73J) became equivalent. In particular the result ([mJ) still holds for the matrices ([73] ) 
up to trivial replacement N — > l±. However, taking into account of even the triple events for the model (fZ3|) , turns 
out to be much more complicated problem than for the model (|5^). 

First of all, the Renormalization Group so successful for the PRBM of the previous section does not work for the 
model (|7^). One may try again to perform the gradual diagonalization (68),(|69|) for the matrices (|73|). However, 

m±(!i-i)' As a result, the orthogonal 



now this procedure should start from the most important sub-diagonal M, t 



transformation ( |68[ ) in this case will change significantly statistics of the remaining part of the matrix. For example, 
if Si — Ei+i ~ Mi.i + /, one has for two elements from the "remaining part" 



M i+ i^ i+ i 
M i+hi 



otiMi+i^+i 
-a 2 Mi + i } i + i 



(74) 



because of Mi + i i <§; Mi + i^ + i due to the (]7^). Suppose for example that ct\ = ot2 = l/V^- This means that due 
to the equation (Q) instead of one large matrix element Mi + i^ + i we have obtained two of them. Each of these two 
matrix elements is \/2 times smaller than the original one, but their sum which in fact is important (see (|63|),(64)) is 
\/2 larger. Thus effectively the secondary interactions in (73) are enhanced compared to those for the PRBM. 
Moreover, the accurate calculation of the triple contribution to the participation ratio for the model (|73|) gives 



P = 1 - V2tt- 



1 



G 



const 



eff 



1 



G 
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G eff = 



G 



Hh) 



(75) 



For short we skip the calculation of the const ~ 1. One sees, that the triple collision contribution compared to the 
PRBM result ( |7C| ) is enhanced by the additional factor ln(G e //). In order to find this logarithm one has to consider 
the effective matrix element of the two stage interaction 



(i\M eff \j) =J2 

k 



M ik M, 



ik ^kj 
£k 



(76) 



The contribution of such matrix element to the participation ratio (see (62 6^)) is proportional to SP ~ — \(i\M e f f\j)\ 
and the averaging over intermediate energy Sk naturally gives the logarithm. As one may show for usual PRBM this 
logarithm does not appear due to the competition between the effective second order matrix element (76) and the 



direct one Mi. 
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The higher order corrections to ( |75| ) may be also shown to have the form const x G^T 1 (In G e //) n . Unfortunately, 
the explicit summation of these logarithmic corrections, even if possible, will hardly teach us anything. On the 
one hand, if G~^ ln(G e //) <C 1 the k pr is also small. On the other hand, the interesting case G~ff ln(G e //) ~ 1 
corresponds to ln(G e //) ~ 1 also and the large logarithm approximation ([75]) does not work in this case. Thus, the 
only conclusion we may draw about the inverse participation ratio and moreover the localization length for the model 
( f73| ) is that they are some unknown functions of G~^ = ln(Zi)/G ([75]). 

Finally, we have the convincing explanation why the chaotic mixing of many states for the model ( |73| ) (as well as 
(f7l|) , (|72]) ) should be enhanced compared to that for the PRBM model (]58|) , (|59|) . However, we have no idea how even 
to estimate analytically the magnitude of this mixing. Therefore, the last possibility we have in hands is to perform 
the numerical simulations. 

The main limitation for our numerical procedure is the impossibility to diagonalizc the very large random matrices. 
Namely, all our data below have been found for matrices of the size L = 100 — 1000 or sometimes L < 2500 (we use 
the notation L for the total size of the matrix in order to stress that it is directly associated with the length of our 
sample). In order not to take into account the boundary effects (for i ~ l\ and L — i ~ Zi, where i is the vector index) 
we have considered the coherent propagation in the periodic ld-sample. The periodic generalization of the model ( |73| ) 
reads 

-G<M nn <G , (77) 

F(i - j) + F(i - j -L) + F{i ~j + L) 



i 3 j MijM kn = {8j k 8 m + dj n d ik ) 
F(x) = 



(h - \x\)~ 2 if < |z| < h 
otherwise 

For comparison of the results for inverse participation ratio U pr with those for more simple model of the previous 
section we have also performed some calculations for the same periodic matrices ( f77| ) but with 

F(x) =x~ 2 , < |at| < h ■ (78) 
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FIG. 1. Participation ratio P for the model (77) - triangles and the simplified power law model (78) - stars. The upper, 
middle and bottom curves corresponds to different strengths of the diagonal G = 10., 5., 2.5 . 
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The participation ratio P averaged over different eigenvectors as well as over many (10 — 50) realizations of the 
random matrix ( f77| , f78|) is shown in the fig. 1 as a function of l\. It is natural to use the logarithmic scale for both P 
and 1%. We have performed the calculations for three values of the strength of the diagonal G = 2.5, 5., 10.. In order 
to get rid of the edge effects close to the border of the energy zone we have taken into account only the vectors with 
eigenvalue < 0.8G for G = 5., 10. and |£j| < 0.7G for G — 2.5. The first point h = 1 coincides for both models as 
may be seen from (|77| , [78| ). However we are interested in the case l\ S> 1. 

The few comments on the figure 1 are of order. 

• The log-log curves for the model ((78|) (solid lines) shows no deviation from the straight line in agreement with 
©■ 

• The slope of both curves for the realistic model ( J7T| ) and for the toy one ([78]) grows with G. This is again in 
agreement with ( ffp ) and with what one may expect from, e.g., ([75]), but disagrees with what one should expect 
for the model of W. 

• The participation ratio for realistic model ( f77| ) decays for large Zi much faster than those for (|78| ) in accordance 
with our wave- handing discussion (75 7^) of the role of satellite peaks in ([72]). 



• The dependence of ln(P) on ln(Zi) for the model (|77|) is not linear at least for l\ < 25 as may be seen from 
the figure. We have fitted these data by a parabola (dashed lines). However, our computations are not enough 
to say with sure either say P ~ exp(— const x ln(Zi) 2 ) for large h or asymptotically P ~ i~ const anc j ^he finite 
curvature of the dashed lines in the fig. 1 accounts only for the small l\ effects? 



The calculation of P was carried out straightforwardly in accordance with the definition (^6|,|62|). The calculation 
of the localization length for circular matrices due to the identification of indices i + L = i requires a little more care. 
Let us consider our circular vector index i literally and define the radius of the "center of mass" for the eigenvector 



R- 



x 2 + y 2 

8=1 



' COS 
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^H 2 sm(^ 



(79) 



where on are the normalized eigenvector components. This R naturally reduces to R = 1 for the very localized 
(pointlike) states and to (R) = for the delocalized ones. One may easily see that for very large matrices the 
deviation of (R) from unit is proportional to the mean squared size of the wave packet I 2 — ((« — «o) 2 )- Therefore for 
explicit definition of h we have used the formula 



l2 = y/((i-(i)) 2 ) = ^\f 



Lll-(R) 



1 + 



(80) 



Here the brackets (...) stands for averaging over the many eigenvectors as well as many realizations of the random 
matrix (fni). In practice we have done a few runs of computations for each Zi with different values of the overall size 
of the matrices L in order to fit out the ~ (l 2 /L) 2 finite size corrections. The finite size effects for l ipr (fig. 1) are 
naturally much weaker ~ exp(— const Ljl%). 

The enhancement factors for the coherent propagation Z2/Z1 found for three different values of the strength of the 
diagonal G — 2.5,5., 10. are shown in the figure 2. It is to be noted that G measures the relative strength of the 
interaction (0)-(0): G ~ U/t for weak coupling and G ~ t/U for strong coupling. 
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FIG. 2. The factor of enhancement for the coherent propagation I2/I1 as a function of l\ for few values of G. Both axes given 
in logarithmic scale. The 7 is the slope of the curve at the asymptotics £1 —* 00 (81). 

As one can see from the fig. 2 the deviations from linear behavior in the logarithm-logarithm plot for the I2/I1 
turn out to be sufficiently weaker than for the participation ratio P. The solid curves on the figure corresponds to 
the least squares fitting with the formula 

y = a + 73; + cexp(— x) , (81) 
y = ln(Z 2 /Zi) , x = ln(Zi) . 



This approximation corresponds to 



h/h-lj [l + f-) . (82) 
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= 10. 
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= 0.56 ±0.04 
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= 5. 
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= 0.76 ± 0.04 
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= 2.5 


7 


= 0.95 ±0.04 



Thus, we have conjectured that the preassymptotic corrections decrease like Z 1 . The values of the exponent 7 for 
different magnitude of the diagonal G are 

1= 10. 7 = 0.56 ±0.04 , 

(83) 

In fact these three 7 may be considered as a main result of our paper. The pure statistical errors for this 7-s are of 
the order of a few percent and may be further improved easily. The main physical problem however is the proper 
choice of the fitting function ( |gl| ) . For example the data presented in the fig. 2 still allows one to use y = a + 72; + cx 2 
instead of (jsil). 

First of all, from the result (|3|) one may conclude with sure that the exponent 7 depends on the strength of the 
diagonal G. However, with only three points in hand we were afraid to fit the dependence 7(G) by some smooth 
function. The physical conditions, evidently consistent with the data (p3]), for the function 7(G) are 7(G = 00) = 
and 7(G = 0) = const ~ 1. Although the case G « 1 does not correspond to any mapping of the original TIP 
problem as we have said after the eq. ( |60| ) . Also the connection between G, and U and t is described in eq. © and 
below. 
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The only rigorous way to confirm numerically the equation (|82j) and to exclude the more exotic dependence of I2 
on li is to consider larger l\. On the other hand, in using of the formula (|8(]) we assume that the total size of the 
sample is large compared to the coherent propagation length L » t Therefore below we would like to develop the 
method which in principle allows one to consider the coherent propagation with only the small matrices L ~ l 2 or 
even L -C h in hand. 

Let us define the new function 



F(L, h) = (y/1 - R(L, h) ) = F{L/l 2 ) . (84) 

Here the last equality will be our scaling hypothesis. Due to ((7^), (80) F = 1 for the very small sample L <§; l 2 and 
F ~ I2J L for L 3> h- For L ~ l 2 the deviation of F from trivial F = 1 will measure the typical warp of the wave 
function on the circle. One should expect that this typical inhomogeneity of the wave function as well as the averaged 
shift of center of mass (R(L,li)) will depend simply on the ratio of the size of the sample L and the two-particle 
localization length 1%. Thus, it is natural to suppose that for L^> l\ the function F{x) is universal in the sense that 
all dependence on li in the r.h.s. of ( p4|) is hidden in hih)- One still may have some doubts in this quite natural 
conjecture due to, e.g., the multifractal nature of the wave function. However, we may easily control its validity within 
our numerical calculations. 

One may fix the value of F and solve numerically the equation for L 

F(L c ,h) = F{x c ) = F c , (85) 
L c = x c (F c ) x l 2 (h) . 

Here only the overall normalization factor x c depends on the value of F c . Of course, the scaling behavior (pi[),(|85|) 
should be violated at L ~ l\. However, because we have seen from fig. 2 that for large single particle localization 
length l 2 S> h one may hope that (|85| ) will be still valid even for l\ <C L <C h- Physically, this means that we are 
trying to find the manifestations of the finite coherent propagation length l 2 in the sample much smaller than this 
length. 
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The dependence of the \n(L c ) on the ln(7i) for G — 5. and for different values F c = 0.2,0.3,0.4,0.5,0.6,0.7,0.8 is 
presented in the hg. 3. At least six of this seven curves looks quite parallel which may be considered as the confirmation 
of the scaling hypothesis (|84|) . Also it is natural that the corrections to scaling as it may be seen from the figure are 
stronger for larger values of F c . In order to take into account at least the corrections to scaling proportional to ~ 
we like in (|8l]) have fitted our results by 

Vi = ^ + (1 + *y)x + Ci exp(-x) , (86) 
Vi=ML {F c )) , x = ln(h) , 

where ai and c$ differs for different F c but 7 is the same for all seven curves of fig. 3. Finally the joint fitting ( |8^ ) 
gives 

7 = 0.83 ±0.04 for G = 5. , (87) 

in agreement with (|8^). The main advantage of this (|S4"|)-([S6|) indirect method of determination of 7 is that we have 
reached the value l\ = 90 (I2 ~ 1300) by considering the L x L matrices with L < 1500. For the direct method 
(|S0|)-(|8"2|) for G = 5. in order to reach even l\ = 25 we have to invoke the much larger matrices with L ~ 2500. 

Of course, due to the logarithmic scale the progress made by going from fig. 2 to fig. 3 may look not so impressive 
and the result (|87f ) due to the scaling hypothesis ( jS4| ) may be considered as model dependent. Nevertheless we 
conclude from (|87f ) that we see no evidence for the violation of the simple power law behavior (^||) for 1%, 



VII. CONCLUSIONS 



In the present paper we have tried to revise the issue of the enhancement of coherent propagation of two interacting 
particles in random potential predicted by D.L. Shepelyansky Q . 

First of all, we definitely see the enhancement but it is of the different nature and has the different functional 
dependence from those usually considered. The existence of the enhancement itself is proved by considering the 
analytical estimate of the matrix element of interaction potential (p3|),(53),(p5|). Effectively the particle-particle 
interaction turns out to be enhanced in logarithm of the Anderson localization length times ln(Zi). Moreover, after 
summation of high order contributions in this effective interaction this logarithm is most likely exponentiated so that 
h/h ~ exp{7ln(/i)} , where 7 = ^(U/t) . 

Unfortunately the functional dependence of the two particle localization length (jS2|) was found only numerically. 
Therefore we still can not completely exclude some other, rather exotic, dependence of I2 on h. The equation j82| ) 
finds also some support in the very similar behavior of the inverse participation ratio li pr ( |70|) which was found 
analytically for the simplified model with the short sample. 

The results ( |82|) , ( |70|) were obtained by mapping of the original two-particle problem onto some random matrix 
model (p58|) , (|59"|) , (|82|) , ([73]) . This means that we are able to explain only qualitatively the behavior of £2- Therefore it is 
quite desirable if somebody will found the same behavior of the coherent propagation length as in our equation ( [82] ) in 
the direct calculation with the Hamiltonian (^) (some numerical results supporting ( |S2"| ) may be found in To this 
end the numerical method of ref. || seems to be very useful. The method of investigation of the effect of localization 
in the samples not large compared to the localization length described at the end of Section 6 being combined with 
the numerical method of j^] may also allow one to consider the larger range of variation of the coherent propagation 
length. 

In general, may be the most interesting result of our paper is that we have found the nontrivial structure (or hierarchy 
of the elements) of the matrix of interaction (^3|) , (|53|) , (p5|) in the basis of noninteracting states (|l7j) , (|33|) , (|4l]) . Just 
due to this hierarchy the matrix models we have to consider differ so drastically from those investigated before. 

In this paper we consider only the interacting particles moving in exactly ld-random potential. It will be very 
interesting to generalize our approach for interacting particles in weak two and three dimensional random potential, 
which are currently the subject of intensive investigation |25|. 
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